# JCR REPLICATION for "How Civilian Loyalties Shape Rebel-led Victimization of Rebel Constituencies" #
# ILAYDA B. ONDER #
# APPENDIX ANALYSIS #
# MAY 27, 2024 #
# This .R file generates the tables and figures in the appendices.

load("Analysis Data.RData") # load analysis data

# Appendix 2: Mccrary test----
#install.packages("rdd")

library(rdd)
DCdensity(df_analysis$margin_kurd, cutpoint = 0, plot = TRUE)
title(xlab = 'Incumbent Party Vote Margin', ylab = 'Density', main = "McCrary Test")
mtext("p = 0.09", side = 3)

# Appendix 3: Table A1 > Table 4 with control variale coefficients----
#install.packages("rdd")
#install.packages("rddtools")
#install.packages("rdrobust")

library(rdd) # package for RDD analysis
library(rddtools) # package for RDD analysis
library(rdrobust) # package for RDD analysis

covs1 = "govrep_vil_log" # control variables for simpler models
covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + 
pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + 
y2014 + y2015 + y2016 + y2017 + y2018" # control variables for extended models

publicworker <- rdd_data(y = publicworker, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data
collaborator <- rdd_data(y = collaborator, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data
politician <- rdd_data(y = politician, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0) # create RDD data

# Run models
f3 <- rdd_reg_lm(rdd_object = publicworker, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 3 of Table 4
se3 <- sqrt(diag(vcovHC(f3, type = "HC1"))) # store robust standard errors

f6 <- rdd_reg_lm(rdd_object = collaborator, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 6 of Table 4
se6 <- sqrt(diag(vcovHC(f6, type = "HC1"))) # store robust standard errors

f9 <- rdd_reg_lm(rdd_object = politician, slope = "separate", order = 1, bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs2, covar.opt = list(strategy = "include")) # Model 9 of Table 4
se9 <- sqrt(diag(vcovHC(f9, type = "HC1"))) # store robust standard errors

# Export Table A.1.
stargazer::stargazer(f3, f6, f9, 
                     type = "text",
                     se = list(se3, se6, se9), 
                     omit = c("x", "x_right", "y2014", "y2015"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", 
                                          "Violent Government Repression (number of incidents)",
                                          "Non-Violent Government Repression (number of incidents)",
                                          "Frequency of armed clashes prior to 2012", 
                                          "Casualties from armed clashes prior to 2012", 
                                          "Extrajudicial killings/political assassinations prior to 2004", 
                                          "Insurgent casualties prior to 2012", 
                                          "Number of voters in 2014 Municipal Elections", 
                                          "Insurgent recruits prior to 2012", 
                                          "District population", 
                                          "Urban population", 
                                          "Kurdish political party vote share in elections prior to 2014", 
                                          "% of Kurdish population", 
                                          "Distance to capital city", 
                                          "Altitude of district centers", 
                                          "Rebellion in the 1920s", 
                                          "Literacy rate", 
                                          "Level of socioeconomic development",
                                          "Year 2016",
                                          "Year 2017",
                                          "Year 2018",
                                          "Constant"),
                     column.labels = c("Public Workers", "Traitors", "Pro-Gov Politicians"),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     title = "Table 4 with Control Variable Coefficients",
                     no.space = T)


# Appendix 4: Table A2 > Calonico-Cattaneo-Titiunik (2014) bandwidth selection----
f1 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(publicworker)$bws[1], 
                 covar.opt = list(strategy = "include"))
se1 <- sqrt(diag(vcovHC(f1, type = "HC1")))

f2 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(publicworker)$bws[1], 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se2 <- sqrt(diag(vcovHC(f2, type = "HC1")))

f3 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(publicworker)$bws[1], 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se3 <- sqrt(diag(vcovHC(f3, type = "HC1")))

f4 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(collaborator)$bws[1], 
                 covar.opt = list(strategy = "include"))
se4 <- sqrt(diag(vcovHC(f4, type = "HC1")))

f5 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(collaborator)$bws[1], 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se5 <- sqrt(diag(vcovHC(f5, type = "HC1")))

f6 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(collaborator)$bws[1], 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se6 <- sqrt(diag(vcovHC(f6, type = "HC1")))

f7 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(politician)$bws[1], 
                 covar.opt = list(strategy = "include"))
se7 <- sqrt(diag(vcovHC(f7, type = "HC1")))

f8 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(politician)$bws[1], 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se8 <- sqrt(diag(vcovHC(f8, type = "HC1")))

f9 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_cct_estim(politician)$bws[1], 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se9 <- sqrt(diag(vcovHC(f9, type = "HC1")))

stargazer::stargazer(f1, f2, f3, f4, f5, f6, f7, f8, f9, 
                     se = list(se1, se2, se3, se4, se5, se6, se7, se8, se9), 
                     type = "text",
                     omit = c("x", "x_right", "y2014", "y2015"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", 
                                          "Violent Government Repression (number of incidents)",
                                          "Non-Violent Government Repression (number of incidents)",
                                          "Frequency of armed clashes prior to 2012", 
                                          "Casualties from armed clashes prior to 2012", 
                                          "Extrajudicial killings/political assassinations prior to 2004", 
                                          "Insurgent casualties prior to 2012", 
                                          "Number of voters in 2014 Municipal Elections", 
                                          "Insurgent recruits prior to 2012", 
                                          "District population", 
                                          "Urban population", 
                                          "Kurdish political party vote share in elections prior to 2014", 
                                          "% of Kurdish population", 
                                          "Distance to capital city", 
                                          "Altitude of district centers", 
                                          "Rebellion in the 1920s", 
                                          "Literacy rate", 
                                          "Level of socioeconomic development",
                                          "Year 2016",
                                          "Year 2017",
                                          "Year 2018",
                                          "Constant"),
                     column.labels = c(rep("Public Workers", 3), rep("Traitors", 3), rep("Pro-Gov Politicians", 3)),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     no.space = T,
                     title = "Alternative Optimal Bandwidth Calculation")



# Appendix 5: Sensitivity to bandwidth selection ----
covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + y2014 + y2015 + y2016 + y2017 + y2018"
bandwidths <- seq(from = 0.05, to = 0.2, by = 0.005)

publicworker <- 
  rdd_data(y = publicworker, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0)
collaborator <- 
  rdd_data(y = collaborator, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0)
politician <- 
  rdd_data(y = politician, x = margin_kurd, data = df_analysis, covar = df_analysis, cutpoint = 0)

## Setup some variables that can be used to store the estimates and confidence intervals
est1 <- rep(NA, length(bandwidths))
ser1 <- rep(NA, length(bandwidths))
est2 <- rep(NA, length(bandwidths))
ser2 <- rep(NA, length(bandwidths))
est3 <- rep(NA, length(bandwidths))
ser3 <- rep(NA, length(bandwidths))

for (i in 1:length(bandwidths)){
  ## Estimate the RD model with the appropriate bandwidth, for instance)
  rd_est_placebo <- rdd_reg_lm(rdd_object = publicworker, bw = bandwidths[i], slope = "separate", order = 1, 
                               covariates = covs2, 
                               covar.opt = list(strategy = "include"))
  rd_est_placebo2 <- rdd_reg_lm(rdd_object = collaborator, bw = bandwidths[i], slope = "separate", order = 1,
                                covariates = covs2, 
                                covar.opt = list(strategy = "include"))
  rd_est_placebo3 <- rdd_reg_lm(rdd_object = politician, bw = bandwidths[i], slope = "separate", order = 1,
                                covariates = covs2, 
                                covar.opt = list(strategy = "include"))
  
  est1[i] <- rd_est_placebo$coefficients[2]
  ser1[i] <- sqrt(diag(vcovHC(rd_est_placebo, type = "HC1")))[2]
  
  est2[i] <- rd_est_placebo2$coefficients[2]
  ser2[i] <- sqrt(diag(vcovHC(rd_est_placebo2, type = "HC1")))[2]
  
  est3[i] <- rd_est_placebo3$coefficients[2]
  ser3[i] <- sqrt(diag(vcovHC(rd_est_placebo3, type = "HC1")))[2]
  
  est = c(est1, est2, est3)
  ser = c(ser1, ser2, ser3)
}

## create a data.frame object with, as variables, cut, est and se
model.name = c(rep("Public Workers", 31), rep("Traitors", 31), rep("Pro-Government Politicians", 31))
out <- data.frame(est,ser,bandwidths, model.name)

## calculate the confidence intervals
out$lo <- out$est - 1.96*out$se
out$hi <- out$est + 1.96*out$se

#install.packages("ggthemes")
#install.packages("lemon")
#install.packages("ggplot2")
#install.packages("tidyr")
#install.packages("dplyr")
#install.packages("tidyverse")
#install.packages("ggthemr")

library(ggthemes)
library(lemon)
library(ggplot2)
library(tidyr)
library(dplyr)
library(tidyverse)
library(ggthemr)
ggthemr::ggthemr("fresh")
out %>%
  mutate(across(model.name, factor, levels = c("Public Workers", "Traitors", "Pro-Government Politicians"))) %>%
  ggplot(aes(x= bandwidths, y= est)) +
  facet_wrap(~ model.name, scales = "free") +
  geom_point(size = 1) +
  geom_linerange(aes(ymin=lo,ymax=hi), size = 1) +
  geom_hline(yintercept = 0, size = 1) +
  scale_color_manual(values = c('black','red')) +
  xlab("Bandwidth") +
  ylab("Estimated Coefficient on Treatment Variable") +
  theme(plot.background = element_rect(color=NA),
        panel.grid.major.y = element_blank(),
        legend.position = "none",
        axis.ticks.length = unit(1,"mm"),
        strip.text.x = element_text(size = 18),
        axis.title = element_text(size = 18),
        legend.title = element_text(size = 18),
        axis.text = element_text(size = 18))



# Appendix 6: Sensitivity to polynomial order----

covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + y2014 + y2015 + y2016 + y2017 + y2018"
poli <- seq(from = 1, to = 4, by = 1)

## Setup some variables that can be used to store the estimates and confidence intervals
est1 <- rep(NA, length(poli))
ser1 <- rep(NA, length(poli))
est2 <- rep(NA, length(poli))
ser2 <- rep(NA, length(poli))
est3 <- rep(NA, length(poli))
ser3 <- rep(NA, length(poli))

for (i in 1:length(poli)){
  ## Estimate the RD model with the appropriate cut-off here (using cutpoint[i], for instance)
  rd_est_placebo <- rdd_reg_lm(rdd_object = publicworker, bw = rdd_bw_ik(kernel = "Normal", publicworker), slope = "separate", covariates = covs2, covar.opt = list(strategy = "include"), order = poli[i])
  rd_est_placebo2 <- rdd_reg_lm(rdd_object = collaborator, bw = rdd_bw_ik(kernel = "Normal", collaborator), slope = "separate", covariates = covs2, covar.opt = list(strategy = "include"), order = poli[i])
  rd_est_placebo3 <- rdd_reg_lm(rdd_object = politician, bw = rdd_bw_ik(kernel = "Normal", politician), slope = "separate", covariates = covs2, covar.opt = list(strategy = "include"), order = poli[i])
  
  est1[i] <- rd_est_placebo$coefficients[2]
  ser1[i] <- sqrt(diag(vcovHC(rd_est_placebo, type = "HC1")))[2]
  
  est2[i] <- rd_est_placebo2$coefficients[2]
  ser2[i] <- sqrt(diag(vcovHC(rd_est_placebo2, type = "HC1")))[2]
  
  est3[i] <- rd_est_placebo3$coefficients[2]
  ser3[i] <- sqrt(diag(vcovHC(rd_est_placebo3, type = "HC1")))[2]
  
  est = c(est1, est2, est3)
  ser = c(ser1, ser2, ser3)
}

## create a data.frame object with, as variables, cut, est and se
model.name = c(rep("Public Workers", 4), rep("Traitors", 4), rep("Pro-Government Politicians", 4))
out <- data.frame(est,ser,poli, model.name)
out$truepoli <- out$poli==1

## calculate the confidence intervals
out$lo <- out$est - 1.96*out$se
out$hi <- out$est + 1.96*out$se

## plot
library(ggthemes)
library(lemon)
library(ggplot2)
library(tidyr)
library(dplyr)
library(tidyverse)
library(ggthemr)
ggthemr::ggthemr("fresh")
out %>%
  mutate(across(model.name, factor, levels = c("Public Workers", "Traitors", "Pro-Government Politicians"))) %>%
  ggplot(aes(x= poli, y= est,color=truepoli)) +
  facet_wrap(~ model.name, scales = "free") +
  geom_point(size = 1) +
  geom_linerange(aes(ymin=lo,ymax=hi), size = 1) +
  geom_hline(yintercept = 0, size = 1) +
  scale_color_manual(values = c('black','red')) +
  xlab("Polynomial Order") +
  ylab("Estimated Coefficient on Treatment Variable") +
  theme(plot.background = element_rect(color=NA),
        panel.grid.major.y = element_blank(),
        legend.position = "none",
        axis.ticks.length = unit(1,"mm"),
        strip.text.x = element_text(size = 18),
        axis.title = element_text(size = 18),
        legend.title = element_text(size = 18),
        axis.text = element_text(size = 18))



# Appendix 7: Models with only incidents included in existing datasets (UCDP, ACLED, GED)----
load("Only GTD ACLED UCDP.RData") # load data containing existing datasets

library(rdd)
library(rddtools)
library(rdrobust)

covs1 = "govrep_vil_log" 
covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + y2014 + y2015 + y2016 + y2017 + y2018"

publicworker <- 
  rdd_data(y = publicworker, x = margin_kurd, data = df_existing, covar = df_existing, cutpoint = 0)
collaborator <- 
  rdd_data(y = collaborator, x = margin_kurd, data = df_existing, covar = df_existing, cutpoint = 0)
politician <- 
  rdd_data(y = politician, x = margin_kurd, data = df_existing, covar = df_existing, cutpoint = 0)
bystander <- 
  rdd_data(y = bystander, x = margin_kurd, data = df_existing, covar = df_existing, cutpoint = 0)

f1 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covar.opt = list(strategy = "include"))
se1 <- sqrt(diag(vcovHC(f1, type = "HC1")))

f2 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se2 <- sqrt(diag(vcovHC(f2, type = "HC1")))

f3 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se3 <- sqrt(diag(vcovHC(f3, type = "HC1")))

f7 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covar.opt = list(strategy = "include"))
se7 <- sqrt(diag(vcovHC(f7, type = "HC1")))

f8 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se8 <- sqrt(diag(vcovHC(f8, type = "HC1")))

f9 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se9 <- sqrt(diag(vcovHC(f9, type = "HC1")))

stargazer::stargazer(f1, f2, f3, f7, f8, f9, 
                     type = "text",
                     se = list(se1, se2, se3, se7, se8, se9), 
                     keep = c("D", "govrep_vil_log", "govrep_non_log"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", "Government Violence", "Government Repression"),
                     column.labels = c(rep("Public Workers", 3), rep("Pro-Gov Politicians", 3)),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     title = "Models with Incidents Included in the UCDP-GED, GTD and ACLED")

# Appendix 10: Additional Analyses Addressing the Potential Endogeneity Issues ----

# Table A.4
load("Tribal Populations Data Appendix.RData") # load additional data for analysis
df_tribalanalysis = merge(df_analysis, appendix_tribal, by = c("province", "district"), all.x = T)

library(rdd)
library(rddtools)
library(rdrobust)

covs1 = "govrep_vil_log" 
covs2 = "govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + y2014 + y2015 + y2016 + y2017 + y2018"

publicworker <- 
  rdd_data(y = publicworker, x = tribe_uprising_per, data = df_tribalanalysis, covar = df_tribalanalysis, cutpoint = 50)
collaborator <- 
  rdd_data(y = collaborator, x = tribe_uprising_per, data = df_tribalanalysis, covar = df_tribalanalysis, cutpoint = 50)
politician <- 
  rdd_data(y = politician, x = tribe_uprising_per, data = df_tribalanalysis, covar = df_tribalanalysis, cutpoint = 50)

f1 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covar.opt = list(strategy = "include"))
se1 <- sqrt(diag(vcovHC(f1, type = "HC1")))

f2 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se2 <- sqrt(diag(vcovHC(f2, type = "HC1")))

f3 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se3 <- sqrt(diag(vcovHC(f3, type = "HC1")))

f4 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covar.opt = list(strategy = "include"))
se4 <- sqrt(diag(vcovHC(f4, type = "HC1")))

f5 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se5 <- sqrt(diag(vcovHC(f5, type = "HC1")))

f6 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se6 <- sqrt(diag(vcovHC(f6, type = "HC1")))

f7 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covar.opt = list(strategy = "include"))
se7 <- sqrt(diag(vcovHC(f7, type = "HC1")))

f8 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs1, 
                 covar.opt = list(strategy = "include"))
se8 <- sqrt(diag(vcovHC(f8, type = "HC1")))

f9 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se9 <- sqrt(diag(vcovHC(f9, type = "HC1")))


stargazer::stargazer(f1, f2, f3, f4, f5, f6, f7, f8, f9, # generate Table A.4.
                     type = "text",
                     se = list(se1, se2, se3, se4, se5, se6, se7, se8, se9), 
                     keep = c("D", "govrep_vil_log", "govrep_non_log"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", "Government Violence", "Government Repression"),
                     column.labels = c(rep("Public Workers", 3), rep("Traitors", 3), rep("Pro-Gov Politicians", 3)),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     title = "Models with Dissident Tribal Populations as the Treatment Variable")


# Table A.5
load("Incumbency Data Appendix.RData") # load additional data for analysis
df_incumbencyanalysis <- merge(df_analysis, appendix_incumbency, by = c("province", "district"), all.x = TRUE)
summary(m1 <- lm(publicworker ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log, data = df_incumbencyanalysis))
te1 <- sqrt(diag(vcovHC(m1, type = "HC1")))

summary(m2 <- lm(publicworker ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco,
                 data = df_incumbencyanalysis))
te2 <- sqrt(diag(vcovHC(m2, type = "HC1")))

summary(m3 <- lm(collaborator ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log, data = df_incumbencyanalysis))
te3 <- sqrt(diag(vcovHC(m3, type = "HC1")))

summary(m4 <- lm(collaborator ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco,
                 data = df_incumbencyanalysis))
te4 <- sqrt(diag(vcovHC(m4, type = "HC1")))

summary(m5 <- lm(politician ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log, data = df_incumbencyanalysis))
te5 <- sqrt(diag(vcovHC(m5, type = "HC1")))

summary(m6 <- lm(politician ~ akpcandidate2014winner2009 + winner2009run2014 + govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco, 
                 data = df_incumbencyanalysis))
te6 <- sqrt(diag(vcovHC(m6, type = "HC1")))

stargazer::stargazer(m1, m2, m3, m4, m5, m6,
                     type = "text",
                     se = list(te1, te2, te3, te4, te5, te6), 
                     covariate.labels = c("Incumbent Candidate runs for the Incumbent Party", 
                                          "Incumbent Candidate runs in the election", 
                                          "Violent Government Repression (# of incidents)",
                                          "Non-Violent Government Repression (# of incidents)",
                                          "Frequency of armed clashes prior to 2012", 
                                          "Casualties from armed clashes prior to 2012", 
                                          "Extrajudicial killings/political assassinations prior to 2004", 
                                          "Insurgent casualties prior to 2012", 
                                          "Number of voters in 2014 Municipal Elections", 
                                          "Insurgent recruits prior to 2012", 
                                          "District population", 
                                          "Urban population", 
                                          "Kurdish political party vote share in elections prior to 2014", 
                                          "% of Kurdish population", 
                                          "Distance to capital city", 
                                          "Altitude of district centers", 
                                          "Rebellion in the 1920s", 
                                          "Literacy rate", 
                                          "Level of socioeconomic development",
                                          "Constant"),
                     column.labels = c(rep("Public Workers", 2), rep("Traitors", 2), rep("Pro-Gov Politicians", 2)),
                     dep.var.labels.include = F,
                     #omit.stat = c("f", "ser"),
                     title = "Models with Incumbency Status of the Candidates in 2014 Elections")

# Appendix 11: Controlling for Other Treatments Happening Simultaneously ----
load("Trustee Data Appendix.RData") # load additional data for analysis
df_trusteeanalysis <- merge(df_analysis, appendix_trustee, by = c("province", "district"), all.x = TRUE)
library(rdd)
library(rddtools)
library(rdrobust)

covs2 = "kayyum2016 + govrep_vil_log + govrep_non_log + frequencyclash + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco + y2014 + y2015 + y2016 + y2017 + y2018"

publicworker <- 
  rdd_data(y = publicworker, x = margin_kurd, data = df_trusteeanalysis, covar = df_trusteeanalysis, cutpoint = 0)
collaborator <- 
  rdd_data(y = collaborator, x = margin_kurd, data = df_trusteeanalysis, covar = df_trusteeanalysis, cutpoint = 0)
politician <- 
  rdd_data(y = politician, x = margin_kurd, data = df_trusteeanalysis, covar = df_trusteeanalysis, cutpoint = 0)


f3 <- rdd_reg_lm(rdd_object = publicworker, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", publicworker), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se3 <- sqrt(diag(vcovHC(f3, type = "HC1")))


f6 <- rdd_reg_lm(rdd_object = collaborator, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", collaborator), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se6 <- sqrt(diag(vcovHC(f6, type = "HC1")))


f9 <- rdd_reg_lm(rdd_object = politician, 
                 slope = "separate", 
                 order = 1, 
                 bw = rdd_bw_ik(kernel = "Normal", politician), 
                 covariates = covs2, 
                 covar.opt = list(strategy = "include"))
se9 <- sqrt(diag(vcovHC(f9, type = "HC1")))


stargazer::stargazer(f3, f6, f9, 
                     type = "text",
                     se = list(se3, se6, se9), 
                     keep = c("D", "kayyum2016"),
                     covariate.labels = c("Treatment (Incumbent Party Victory)", "Trustee Appointment"),
                     column.labels = c(rep("Public Workers", 1), rep("Traitors", 1), rep("Pro-Gov Politicians", 1)),
                     dep.var.labels.include = F,
                     omit.stat = c("f", "ser"),
                     title = "Models Controlling for the Assignment of Custodians to Municipalities Won by HDP")


# Appendix 12: Multinomial Logit Models----
load("Multinomial Logit Data Appendix.RData")
#install.packages("nnet")
library(nnet)

multinomial$akp_win <- as.factor(multinomial$akp_win)
summary(mul1 <- multinom(outcome ~ 1 + akp_win + govrep_vil_log + govrep_non_log
                         , 
                         data = multinomial))

summary(mul2 <- multinom(outcome ~ 1 + akp_win + govrep_vil_log + govrep_non_log +
                           frequencyclash  + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco
                         , 
                         data = multinomial))


summary(mul3 <- multinom(outcome ~ 1 + akp_win + govrep_vil_log + govrep_non_log +
                           frequencyclash  + intensityclash + murders + deaths + voters + recruitswe + pop + urbanpop + kurdvote + kurdpop + distance + altitude + rebellion + literacy + soceco +
                           y2016 + y2017 + y2018
                         , 
                         data = multinomial))
install.packages("stargazer")
library(stargazer)
stargazer(mul1, mul2, mul3, # generate Table A.7
          type = "text",
          no.space = T,
          covariate.labels = c("Treatment (Incumbent Party Victory)", 
                               "Violent Government Repression (number of incidents)",
                               "Non-Violent Government Repression (number of incidents)",
                               "Frequency of armed clashes prior to 2012", 
                               "Casualties from armed clashes prior to 2012", 
                               "Extrajudicial killings/political assassinations prior to 2004", 
                               "Insurgent casualties prior to 2012", 
                               "Number of voters in 2014 Municipal Elections", 
                               "Insurgent recruits prior to 2012", 
                               "District population", 
                               "Urban population", 
                               "Kurdish political party vote share in elections prior to 2014", 
                               "% of Kurdish population", 
                               "Distance to capital city", 
                               "Altitude of district centers", 
                               "Rebellion in the 1920s", 
                               "Literacy rate", 
                               "Level of socioeconomic development",
                               "Year 2016",
                               "Year 2017",
                               "Year 2018",
                               "Constant"
          ))




# Predictions: Figure A.4
require(effects)


# Constituency
probs1 = Effect("akp_win", mul2)
preds = as.data.frame(probs1$prob)
lower = as.data.frame(probs1$lower.prob)
upper = as.data.frame(probs1$upper.prob)
preds$value = c("Incumbent Party Loss", "Incumbent Party Victory")
lower$value = c("Incumbent Party Loss", "Incumbent Party Victory")
upper$value = c("Incumbent Party Loss", "Incumbent Party Victory")
preds = gather(preds, condition, measurement, prob.politician:prob.collaborator, factor_key=TRUE)
lower = gather(lower, condition, measurement, L.prob.politician:L.prob.collaborator, factor_key=TRUE)
upper = gather(upper, condition, measurement, U.prob.politician:U.prob.collaborator, factor_key=TRUE)
lower$condition = as.character(lower$condition)
lower$condition = substring(lower$condition, 3, nchar(lower$condition))
names(lower)[3] = "lower"
upper$condition = as.character(upper$condition)
upper$condition = substring(upper$condition, 3, nchar(upper$condition))
names(upper)[3] = "upper"
preds = merge(preds, lower, by = c("value", "condition"))
preds = merge(preds, upper, by = c("value", "condition"))
preds$condition = c("Collaborator", "Politician", "Public Worker",
                    "Collaborator", "Politician", "Public Worker")
preds$condition <- factor(preds$condition, levels = c("Collaborator", "Politician", "Public Worker"))
preds$model = "All Groups"

install.packages("facetscales")
library(facetscales)

scales_y <- list(
  `Cooperation` = scale_y_continuous(limits = c(0, 0.5), breaks = seq(0, 0.5, 0.1)),
  `No engagement` = scale_y_continuous(limits = c(0.3, 1), breaks = seq(0.3, 1, 0.1)),
  `Infighting` = scale_y_continuous(limits = c(0, 0.7), breaks = seq(0, 0.7, 0.1))
)
library(ggplot2)
library(facetscales)
r = ggplot(preds, aes(x = value, y = measurement, ymin = lower, ymax = upper, group = 1, colour = condition)) +
  #facet_wrap( ~ condition, ncol = 3, scales = "free") +
  facet_grid_sc(rows = vars(condition)) +
  ggtitle("") +
  theme_bw() +
  geom_point(size = 3) +
  geom_line(size = 1) +
  geom_errorbar(size = 1, width = 0.05) +
  labs(x = "", 
       y = "") +
  theme(plot.title = element_text(size = 15, hjust = 0.5),
        axis.title = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12, angle = 45, vjust = 1, hjust = 1),
        strip.text = element_text(size = 12), 
        legend.position = "none",
        plot.margin = unit(c(1,1,0.5,1), "cm"))
r

#END ----

